CREST Adapter - #807
Conversation
d6a83b4 to
c231a43
Compare
3e88c36 to
7674f5f
Compare
9c25f3d to
51368be
Compare
9be935e to
eab7647
Compare
An impossible spin multiplicity (e.g. an even multiplicity for a species
with an even electron count) was forwarded all the way to the ESS, where
it surfaced as a cryptic non-retryable Gaussian GL301 InputError and
aborted the species after wasting a full run.
Enforce the parity relation
(total_electrons - net_charge) % 2 == (multiplicity - 1) % 2
in ARCSpecies once charge and multiplicity are finalized, raising a clear
SpeciesError that names the label, electron count, net charge, requested
multiplicity, and the nearest valid multiplicities. This is an inviolable
parity relation, so the guard has zero false positives. It is skipped
gracefully when the electron count is not yet known (e.g. a TS or an
xyz-less species with no perceived Molecule).
A fine transition-state optimization that dead-ends on the Gaussian l9999 "Optimization stopped" / "Number of steps exceeded" oscillation (classified as ['MaxOptCycles', 'GL9999']) was previously escalated only through the Hessian/step ladder (maxcycle=200 -> recalcfc=5 -> calcall -> RFO) and then abandoned via 'all_attempted'. opt=(cartesian) - wired until now only to the GL103/InternalCoordinateError path - was never tried, even though the failure is a redundant-internal-coordinate oscillation that Cartesian coordinates sidestep. This lost chemically-sensible TS guesses (e.g. arcbench reaction_16 C=C + HCl, heuristics-xy conformers 0/1). Insert opt=(cartesian) as an early rung in the MaxOptCycles ladder, right after the cheap maxcycle bump and before the expensive Hessian/algorithm escalation, gated tightly to the fine TS opt (is_ts and fine). Ground-state opts (is_ts=False) and coarse TS opts (fine=False) are unchanged, as is the GL103 path. Cartesian is recorded once (bare 'cartesian' marker, consistent with trsh_keyword_cartesian) and folded into the merged opt route so it is carried through the remaining rungs; the ladder still terminates via 'all_attempted' if Cartesian also fails. Also fixes a latent cartesian-loss in mixed GL103+MaxOptCycles histories, where the opt-route rewrite would previously clobber a standalone opt=(cartesian).
| seed = {'xyz': xyz, 'family': rxn.family} | ||
| constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) | ||
| self.assertIsInstance(constraints, dict) | ||
| self.assertTrue({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'} <= set(constraints)) |
Check notice
Code scanning / CodeQL
Imprecise assert Note
Show autofix suggestion
Hide autofix suggestion
Copilot Autofix
AI 9 days ago
Replace the imprecise boolean assertion with a specific unittest assertion tailored to subset comparison.
Best fix in this file: in arc/job/adapters/ts/heuristics_test.py, method test_get_wrapper_constraints_crest, replace:
self.assertTrue(required_keys <= set(constraints))style check
with:self.assertLessEqual(required_keys, set(constraints))
This preserves exact test logic (subset relationship) while improving failure output (you’ll see the compared sets). No new imports or helper methods are needed.
| @@ -2391,7 +2391,7 @@ | ||
| seed = {'xyz': xyz, 'family': rxn.family} | ||
| constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed) | ||
| self.assertIsInstance(constraints, dict) | ||
| self.assertTrue({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'} <= set(constraints)) | ||
| self.assertLessEqual({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'}, set(constraints)) | ||
| self.assertEqual( | ||
| (constraints['A'], constraints['H'], constraints['B']), | ||
| constraints['angle_atoms'], |
…ny cores On small double-radical PESs (e.g. OH + OH -> H2O + O) the CREST H-abstraction $constrain block pinned only the A-H and H-B distances plus a weak positional restraint, leaving the A...B heavy-heavy separation and the A-H-B angle free. GFN2-xTB metadynamics then collapsed the linear A...H...B seed into a bent A-B minimum, which the post-CREST angle guard correctly rejected. Part 1: for the three-center H-abstraction path, also emit 'distance: A, B, auto' and 'angle: A, H, B, auto' in the $constrain block (xtb supports the angle keyword; see xtb_adapter.py). The XY four-center path (no angle_atoms) is left unchanged. Part 2: skip the CREST guess when the reactive core spans essentially the whole molecule (H_Abstraction with <=1 spectator atom, i.e. <=4 atoms total), where CREST has no conformational DOF to sample; defer to the other TS-search methods.
An l202 'OptOrientation' error (standard orientation / point group changed mid-optimization) is a symmetry glitch that nosymm cures - the same remedy l101/l103 already use. It was mis-grouped with genuine dead-ends in GAUSSIAN_NON_RETRYABLE_KEYWORDS, which short-circuited trsh_ess_job into a non-retryable refusal before the Gaussian remedy branch was reached. Remove 'OptOrientation' from GAUSSIAN_NON_RETRYABLE_KEYWORDS and add 'NoSymm' to the l202 keyword list in determine_ess_status, so trsh_keyword_nosymm fires and injects nosymm exactly once. If nosymm also fails the job terminates via all_attempted (no infinite retry). l101/l103 and other non-retryable keywords are unchanged.
For a reaction whose family was discovered in the reverse direction, the
first-orientation map_rxn attempt fails on get_template_product_order (or on a
missing template map) and returns None, after which map_reaction routinely
recovers via its existing flip fallback and returns a valid map. The two
intermediate logger.error calls in map_rxn therefore fired as false alarms on
reactions that mapped successfully.
Downgrade both to logger.debug. Control flow is unchanged (both branches still
return None as before), so the returned atom map is identical for every
reaction. A genuine total-mapping failure (both orientations exhausted) is
still surfaced at error level by ARCReaction.atom_map ("could not be atom
mapped").
Add a focused test that drives the first-orientation failure via map_rxn and
asserts no ERROR is logged, while map_reaction still returns a valid map.
ARCSpecies.set_dihedral() correctly refuses to rotate a dihedral that spans a linear segment (e.g. a collinear O=C=C cumulene) and returns None. This is a benign no-op, but it was logged at WARNING level and the mapping engine's backbone-alignment loop re-attempted the same linear torsion on every iteration, producing repeated warnings. - set_dihedral(): downgrade the linear-segment message to logger.debug. The is_angle_linear guard and the return None are unchanged. - get_backbone_dihedral_angles(): skip backbone torsions that span a linear segment (new is_torsion_linear helper, mirroring the guard) so they are never fed into the set_dihedral alignment loop. This is a no-op for non-linear torsions and does not change alignment results. Adds focused tests for both behaviors.
Molpro's `memory,Total=N,m` card is allocated PER MPI PROCESS and the job
runs `molpro -n {cpu_cores}`, but set_input_file_memory hard-coded N as
job_memory_gb * 31.25 (== 125/4, an assumed 4 processes) and never divided
by the actual cpu_cores. At higher rank counts the card demanded several
times the granted node memory -> 'Insufficient memory to allocate'.
(a) Compute the card per-process: ceil(job_memory_gb * 125 / cpu_cores),
mirroring Orca/CFOUR/TeraChem. 125 MW == 1 GB, so total memory equals
job_memory_gb for any rank count (previously only correct at 4 ranks).
(b)/(c) Molpro 'Memory' trsh now bounds the requested total by the node
cap (node_mem * job_max_server_node_memory_allocation) instead of
inflating unboundedly, and once at the cap reduces the MPI rank count
(halving cpu_cores, already returned/threaded to set_cpu_and_mem) so the
per-process share grows without exceeding node memory. When it cannot
reduce ranks further it emits a terminal error rather than resubmitting
identically until max_ess_trsh.
Update molpro_test set_input_file_memory expectations (48->37, 8->219,
1->1750 for job_memory_gb=14) and the written memory cards (438->37 at 48
ranks); update/extend trsh_test for the bounded + rank-reducing branch.
The at-cap rank-halving branch set memory = min(memory_gb, max_mem_allocation), which kept the pre-cap total whenever memory_gb was still below the cap. A job that needs ~the full node (e.g. 62 GB, cap 60.8) was then declared dead at a single rank having never actually been offered the cap. Pin the total to the node cap (memory = max_mem_allocation) while reducing ranks so the per-rank share is raised against the full node allocation before giving up. Update trsh_test expectations: the at-cap halving cases now expect memory == max_mem_allocation (the cap) at each step; the focused test's at-cap case uses memory_gb below the cap so the pin lifts the total up to the cap.
…ount) Atomic nitrogen has 7 electrons, so its multiplicity must be even. Ground-state N is 4S (multiplicity=4). The old multiplicity=3 is physically impossible and now raises SpeciesError under the electron-parity guard, erroring the whole molpro_test class. Also update the job_7 expected Molpro input wf,spin cards from spin=2 to spin=3 (Molpro spin = multiplicity - 1).
…d trsh The prior fix divided the memory,Total=N,m card by cpu_cores on the assumption it is per-process. molpro26 on zeus proves it is NODE-TOTAL: a written Total=438 printed 'Total memory per node: 438 MW' while giving only 'Memory per process: 30 MW' (Molpro reserves ~25% for Global Arrays and splits the rest across the -n ranks itself). Dividing by cpu_cores double-divides and starves the node; a 42 GB/12-rank job got 3.5 GB. - molpro.py set_input_file_memory: N = ceil(job_memory_gb * 125), so the node-total card equals the PBS reservation. No cpu_cores division. - trsh.determine_ess_status: the 'A further ... Increase memory to M' branch now captures the target M (last number) as a per-process floor; fixes the always-false 'error not in locals()' guard (-> 'not error'). - trsh.trsh_ess_job molpro Memory branch: when a per-process target is parsed, size the node-total card as ceil(M*1.5*nprocs/0.822)/125 GB so each rank gets >= M MW; keep the memory_gb*3 fallback and the at-cap rank-halving. Fixed stray /128 -> node-total sizing. - Tests: node-total card == ceil(job_memory_gb*125), independent of cpu_cores; added a targeted-retry test for the per-process-target path.
The heuristic TS seed-builder (combine_coordinates_with_redundant_atoms) failed to produce any guess when the H-abstraction acceptor has a linear / cumulene reactive center (e.g. HCCO + THF <=> CH2CO + THF_rad, where the ketene product has a collinear O=C=C backbone). Root cause: in get_new_zmat_2_map, after remove_zmat_atom_0 the reduced zmat_2 of a linear/cumulene fragment carries a collinear (~0 deg) bond angle. zmat_to_xyz cannot place a dihedral that references three collinear atoms (cross product is zero -> divide-by-zero -> NaN), so every atom collapses to the origin and the ARCSpecies build for atom mapping raises a colliding-atoms SpeciesError, which aborted the combine and yielded no TS guesses. That species is used only for connectivity-based atom mapping, so on SpeciesError we rebuild it from a copy of the zmat with collinear angles nudged off exactly 0/180 deg (connectivity preserved). The returned TS coordinates are unaffected: they are built from the original, un-perturbed zmat_2. Non-linear behavior is byte-for-byte unchanged (the fallback only runs when the direct build fails). Adds a regression test exercising a cumulene (O=C=C) acceptor.
…n seed build fails When CREST's heuristic Z-matrix seed construction yields no seed guesses (e.g. a linear/cumulene reactive center such as HCCO in H_Abstraction that the heuristic builder cannot assemble), CREST previously contributed nothing even when another adapter (e.g. AutoTST) had already produced valid TS guesses for the same TS. Add a graceful, opt-in fallback: seed_hub.get_backup_ts_seeds() harvests successful, non-CREST TS guesses already present on rxn.ts_species.ts_guesses (the incore TS adapters run sequentially with CREST last, so earlier guesses are available in-process) and returns them as CREST seeds with empty metadata. The existing get_wrapper_constraints() machinery then re-derives the family reactive-atom / distance constraints from the seed geometry itself, so an externally-built guess is a valid seed without the failed geometry construction. CrestAdapter.execute_incore() only takes this path when its own seed build returns nothing; the normal (seed-succeeds) path is unchanged. A feedback-loop guard (exclude_method='crest') ensures CREST is never seeded from a prior CREST result, and the pre-existing tiny-system / family gates still apply upstream. Adds two focused unit tests to crest_test.py covering seed selection (feedback/failed-guess exclusion, opt_xyz preference, dedup) and the seed-injection -> constraint-derivation -> input-generation path.
| if param.split('_')[0] in ('A', 'AX'): | ||
| if value <= delta: | ||
| new_vars[param] = delta | ||
| elif value >= 180.0 - delta: |
Check warning
Code scanning / CodeQL
Redundant comparison Warning
Show autofix suggestion
Hide autofix suggestion
Copilot Autofix
AI 8 days ago
General fix: remove ambiguity in chained range comparisons by computing explicit bounds from validated inputs, then compare against those bounds. This avoids analyzer confusion and makes intent clear.
Best fix here: in arc/job/adapters/ts/heuristics.py, inside _perturb_collinear_zmat_angles, normalize delta to a safe perturbation magnitude in (0, 90) and compute lower_bound/upper_bound once. Then replace direct comparisons against delta / 180.0 - delta with comparisons against these bounds. This keeps current functionality for normal calls (delta=5.0) and prevents pathological values from creating overlapping conditions.
Changes needed:
- Edit only
_perturb_collinear_zmat_anglesbody (around lines 729–735). - No new imports, methods, or external dependencies.
| @@ -726,13 +726,19 @@ | ||
| Returns: | ||
| dict: A copy of ``zmat`` with collinear bond angles perturbed. | ||
| """ | ||
| perturb = abs(float(delta)) | ||
| if perturb >= 90.0: | ||
| perturb = 89.999 | ||
| lower_bound = perturb | ||
| upper_bound = 180.0 - perturb | ||
|
|
||
| new_vars = dict(zmat['vars']) | ||
| for param, value in zmat['vars'].items(): | ||
| if param.split('_')[0] in ('A', 'AX'): | ||
| if value <= delta: | ||
| new_vars[param] = delta | ||
| elif value >= 180.0 - delta: | ||
| new_vars[param] = 180.0 - delta | ||
| if value <= lower_bound: | ||
| new_vars[param] = lower_bound | ||
| elif value >= upper_bound: | ||
| new_vars[param] = upper_bound | ||
| return {**zmat, 'vars': new_vars} | ||
|
|
||
|
|
… N ranks, was Y GB across M)
…ocks determine_ess_status() mined the Orca SCF block for an initial/final energy pair to run its divergence heuristic. Two ways that crashed on a normally terminated run: 1. The scan for the first SCF iteration energy assumed the line immediately following the incremental-Fock marker holds one. Orca may print a SOSCF header and a separator in between, and the separator holds a single token, so `forward_lines[j + 1].split()[1]` raised IndexError. The loop was also unbounded and could only terminate by raising. 2. `scf_energy_last_iteration` was assigned only inside the 'TOTAL SCF ENERGY' branch but read whenever an initial energy was found, so a log with a parseable iteration line and no 'TOTAL SCF ENERGY' header raised UnboundLocalError. Neither exception is caught by any caller, so either one aborts the whole run. Cherry-picked from fix_orca_ess_status (b79d011).
determine_ess_status() dispatches on an exact software name. OrcaNEBAdapter passed 'orca_neb', matched no branch, and fell through to the terminal `return '', list(), '', ''`, so a normally terminated NEB run carried an empty status rather than 'done'. Since 35d9bd3 the tsg ingestion gate requires ess_succeeded for .log output, so that empty status silently discarded every Orca NEB TS guess. Measured on the 500-reaction benchmark: orca_neb contributed guesses in 10 of 16 reactions before 2026-07-15 and in 0 of 15 after. Override ess_software to 'orca' on OrcaNEBAdapter, and use it at the remaining determine_ess_status call site that still passed job_adapter. Note that Orca NEB now takes part in Orca's ESS error handling like any other Orca job; in particular a 'Syntax' keyword raises JobError, which Scheduler.end_job turns into a re-run. Cherry-picked from fix_orca_ess_status (a0e1a6a).
test_determine_ess_status_orca_scf_energies_unavailable built its two logs as
inline strings. Replace them with fixtures derived from real Orca output:
arc/testing/neb/neb_res.out truncated at two points that leave the SCF block
incomplete, so both carry genuine Orca formatting (the SOSCF header and the
one-token S-O-S-C-F separator that caused the original IndexError) rather than
a hand-written approximation.
No log already under arc/testing/ exhibits either shape -- the corpus was
collected from runs that failed in interesting ways, and a run that stops
before emitting TOTAL SCF ENERGY looks like infrastructure noise -- so the
fixtures are derived rather than captured.
Each pins one half of the fix:
neb_scf_no_iteration_line.out marker + SOSCF header + column header, no
numeric row. Pre-fix this raised IndexError on
the separator; the bounded scan must instead
run off the end without assigning an energy.
neb_scf_no_total_energy.out a parseable iteration row but no TOTAL SCF
ENERGY block, so the initial energy is
assigned and the final one is not. With the
scan fixed but the guard reverted this still
raises, so it isolates the guard.
Verified red with arc/job/trsh.py reverted to its pre-fix state, green with it
restored. Drops the _write_temp_log helper.
Two changes to the TS-guess completion path. The ingestion gate added for QST2 requires a successful ESS run before a .log output is read. QST2 and Orca NEB both write a file called input.log, so the gate caught NEB as well -- but the two differ: a failed QST2 job has nothing usable in its log, whereas NEB writes its geometry before the post-processing that may then fail. Orca NEB contributed guesses in 10 of 16 benchmark reactions before that gate and none of 15 after. Key the exemption on the adapter rather than on a file extension several unrelated adapters share, so a failed QST2 job stays excluded. troubleshoot_ess already knows about tsg jobs -- it records un-troubleshootable errors onto the TSGuess -- but no call site ever passed it one, so that branch was unreachable and a guess lost to a recoverable failure was dropped with only a warning. Call it for an errored TS-guess job, guarded on times_rerun so a persistent failure is attempted once.
Addition of CREST Adapter that complements the heuristic adapter.
This pull request adds support for the CREST conformer and transition state (TS) search method to the ARC project, along with several related improvements and code cleanups. The most important changes include integrating CREST as a TS search adapter, updating configuration and constants, and enhancing the heuristics TS search logic for better provenance tracking and code clarity.
CREST Integration:
JobEnum(arc/job/adapter.py), included CREST in the list of adapters and RMG family mapping, and registered it as a default incore adapter (arc/job/adapters/common.py,arc/job/adapters/ts/__init__.py). [1] [2] [3] [4]arc/job/adapters/ts/crest_test.py).Makefile). [1] [2]Constants and Configuration:
angstrom_to_bohrconversion constant to both Cython and Python constants modules (arc/constants.pxd,arc/constants.py). [1] [2] [3]Heuristics TS Search Enhancements and Refactoring:
arc/job/adapters/ts/heuristics.py). [1] [2] [3] [4]arc/job/adapters/ts/heuristics.py). [1] [2] [3] [4] [5] [6] [7].